---
title: "R18.LandReform_Tables"
---

# Load Packages
```{r}
library(tidyverse)
library(sandwich)
library(stargazer)
```

# Load Data
```{r}
load("Data/data_final.Rdata")
load("Data/gnr_panel.Rdata")
load("Data/ideo_list.Rdata")
```

# Set Covariates (Including Interactions with Each Election Year)
```{r}
covs_list <- list()

covs_list[[1]] <- c("tx_Hi", "election", "tx_Hi:election", "election:as.factor(DistN)", "election:ln_PopDens_70", "election:Cities_10km", "election:ln_Confs_All", "election:PrRecH_Adm3", "election:Isolados1950", "election:PrWheatAr_1972", "election:LandIneqRat", "election:TerrainSlopeIndex", "election:MedAltitude",  "election:Communist_1975", "election:MinDistPCP_Map", "election:log_gnr_cumul_23_Z")
covs_list[[2]] <- c("tx_Hi", "election", "tx_Hi:election", "election:as.factor(DistN)", "election:ln_PopDens_70", "election:Cities_10km", "election:ln_Confs_All", "election:PrRecH_Adm3", "election:Isolados1950", "election:PrWheatAr_1972", "election:LandIneqRat", "election:TerrainSlopeIndex", "election:MedAltitude",  "election:Communist_1975", "election:MinDistPCP_Map")
```

```{r}
elections <- c("AC75", "AR76", "AR79", "AR80", "AR83", "AR85", "AR87")
```

# Generate First Differences between Every Election Year of Interest (1976-87) and Baseline Year (1975)
```{r}
aggregated_fd <- do.call(rbind, ideo_list[1:7]) %>% left_join(ideo_list[[1]], by = c("DICOFRE_Z", "ideology2")) %>% 
  mutate(FD = votes.x - votes.y) %>%
  select(DICOFRE_Z, ideology2, FD, election.x) %>% 
  rename(election = election.x) %>% 
  pivot_wider(c(DICOFRE_Z, election), names_from = "ideology2", values_from = "FD")
```

# Join vote share first difference and cumulative GNR events up to election year to main data set
```{r}
data_models <- data_final %>% 
  left_join(aggregated_fd) %>% 
  left_join(gnr_panel %>% rename(election = pre_election))
```

# Estimate Low LR (App. Table 7ai)
```{r}

models<-list()
j<-1 # full set of covariates

models[[1]] <- lm(as.formula(paste("Communist", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
                data = data_models %>% filter(exp_cat2 == 1))
models[[2]]  <- lm(as.formula(paste("Center", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
               data = data_models %>% filter(exp_cat2 == 1))
models[[3]]  <- lm(as.formula(paste("Left", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
                data = data_models %>% filter(exp_cat2 == 1))
models[[4]]  <- lm(as.formula(paste("Right", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
                data = data_models %>% filter(exp_cat2 == 1))

# Robust standard errors
robust_se1 <- sqrt(diag(vcovHC(models[[1]], type = "HC2")))
robust_se2 <- sqrt(diag(vcovHC(models[[2]], type = "HC2")))
robust_se3 <- sqrt(diag(vcovHC(models[[3]], type = "HC2")))
robust_se4 <- sqrt(diag(vcovHC(models[[4]], type = "HC2")))

# Stargazer output 
stargazer(models,
          se = list(robust_se1, robust_se2, robust_se3, robust_se4),
          dep.var.labels = c("Communist", "Center-Left", "Far Left", "Right"),
          type = "text",
          keep = c("tx_Hi:election"), 
          covariate.labels = c("Treatment, 1976", "Treatment, 1979", "Treatment, 1980", "Treatment, 1983", "Treatment, 1985", "Treatment, 1987"),
          add.lines = list(c("Observations Per Year", length(models[[1]]$residuals)/7, length(models[[2]]$residuals)/7, length(models[[3]]$residuals)/7, length(models[[4]]$residuals)/7)))

```

# Estimate High LR (App. Table 7aii)
```{r}

models<-list()
j<-1 # full set of covariates

models[[1]] <- lm(as.formula(paste("Communist", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
                data = data_models %>% filter(exp_cat2 == 2))
models[[2]]  <- lm(as.formula(paste("Center", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
               data = data_models %>% filter(exp_cat2 == 2))
models[[3]]  <- lm(as.formula(paste("Left", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
                data = data_models %>% filter(exp_cat2 == 2))
models[[4]]  <- lm(as.formula(paste("Right", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
                data = data_models %>% filter(exp_cat2 == 2))

# Robust standard errors
robust_se1 <- sqrt(diag(vcovHC(models[[1]], type = "HC2")))
robust_se2 <- sqrt(diag(vcovHC(models[[2]], type = "HC2")))
robust_se3 <- sqrt(diag(vcovHC(models[[3]], type = "HC2")))
robust_se4 <- sqrt(diag(vcovHC(models[[4]], type = "HC2")))

# Stargazer output 
stargazer(models,
          se = list(robust_se1, robust_se2, robust_se3, robust_se4),
          dep.var.labels = c("Communist", "Center-Left", "Far Left", "Right"),
          type = "text",
          keep = c("tx_Hi:election"), 
          covariate.labels = c("Treatment, 1976", "Treatment, 1979", "Treatment, 1980", "Treatment, 1983", "Treatment, 1985", "Treatment, 1987"),
          add.lines = list(c("Observations Per Year", length(models[[1]]$residuals)/7, length(models[[2]]$residuals)/7, length(models[[3]]$residuals)/7, length(models[[4]]$residuals)/7)))


```

# Estimate LR None (App. Table 7aiii)
```{r}

models<-list()
j<-2 # drop GNR Cumulative Events covariate, as none of these parished had any land expropriated

data_models2 <- data_models %>% 
  mutate(log_gnr_cumul_23_Z = ifelse(exp_cat2 == 0, NA, log_gnr_cumul_23_Z))

models[[1]] <- lm(as.formula(paste("Communist", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
                data = data_models2 %>% filter(exp_cat2 == 0))
models[[2]]  <- lm(as.formula(paste("Center", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
               data = data_models2 %>% filter(exp_cat2 == 0))
models[[3]]  <- lm(as.formula(paste("Left", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
                data = data_models2 %>% filter(exp_cat2 == 0))
models[[4]]  <- lm(as.formula(paste("Right", paste(c(covs_list[[j]]), collapse= " + "), sep = " ~ ")),
                data = data_models2 %>% filter(exp_cat2 == 0))

# Robust standard errors
robust_se1 <- sqrt(diag(vcovHC(models[[1]], type = "HC2")))
robust_se2 <- sqrt(diag(vcovHC(models[[2]], type = "HC2")))
robust_se3 <- sqrt(diag(vcovHC(models[[3]], type = "HC2")))
robust_se4 <- sqrt(diag(vcovHC(models[[4]], type = "HC2")))

# Stargazer output 
stargazer(models,
          se = list(robust_se1, robust_se2, robust_se3, robust_se4),
          dep.var.labels = c("Communist", "Center-Left", "Far Left", "Right"),
          type = "text",
          keep = c("tx_Hi:election"), 
          covariate.labels = c("Treatment, 1976", "Treatment, 1979", "Treatment, 1980", "Treatment, 1983", "Treatment, 1985", "Treatment, 1987"),
          add.lines = list(c("Observations Per Year", length(models[[1]]$residuals)/7, length(models[[2]]$residuals)/7, length(models[[3]]$residuals)/7, length(models[[4]]$residuals)/7)))


```


